#include "mesh.hpp"
#include "topology.hpp"
#include "signalCodes.hpp"
#include "exceptions.hpp"
#include "util.hpp"
#include "meshFunctions.hpp"
#include "macros.hpp"
#include <algorithm>
#include <cassert>
#include <experimental/memory>
#include <utility>

#include "timer.hpp"


namespace hms
{

#ifndef NDEBUG
Mesh::Mesh( Matrix2Xs vertexCoords ) :
	m_vertices { std::move(vertexCoords) }
{}

Mesh::Mesh( Matrix2Xs vertexCoords, Array2Xi vertexIndices, bool ) :
	Mesh( std::move(vertexCoords) )
{
	m_edges.vertices( std::move(vertexIndices) );
}
#endif


/**
 * @brief Constructor for unstructured meshes.
 * Use this ctor when you have vertex coordinates as well as the vertex indices
 * (at either end) of all edges.
 * 
 * @param vertexCoords 
 * @param edgeVertexIndices 
 */
Mesh::Mesh( Matrix2Xs&& vertexCoords, Array2Xi&& edgeVertexIndices ) :
	Mesh( edgesToMesh(
		std::move(vertexCoords),
		std::move(edgeVertexIndices)
	) )
{}

/**
 * @brief Constructor for unstructured meshes.
 * Use this ctor when you have vertex coordinates as well as cells defined by
 * their number of vertices (=number of edges) and vertex indices.
 * VTK and gmsh store vertex indices for cells, so this will be handy for
 * importers.
 * 
 * @param vertexCoords 
 * @param nDomainCellVertices
 * @param domainCellVertices 
 */
Mesh::Mesh(
	Matrix2Xs&& vertexCoords,
	ArrayX<char>&& nDomainCellVertices,
	const Array4Xi& domainCellVertices
) : Mesh( cellVerticesToMesh(
		std::move(vertexCoords),
		std::move(nDomainCellVertices),
		domainCellVertices
	) ) // this effectively chains to another ctor (see function def for more).
{}

/**
 * @brief Constructor for unstructured meshes.
 * Use this ctor when you have
 * vertex coordinates,
 * vertex indices of edges,
 * domain cells defined by edge indices (which also requires the number of edges
 * for each cell).
 * If available, you can pass in vertex indices for cells. They are necessary 
 * to determine cell areas and centroids, but can be reconstructed if not passed
 * in.
 * Cell centroids are also required to determine the orientation of edge normal
 * vectors.
 * 
 * @param vertexCoords 
 * @param edgeVertexIndices 
 * @param nDomainCellEdges 
 * @param domainCellEdges 
 * @param domainCellVertices 
 */
Mesh::Mesh(
	Matrix2Xs&& vertexCoords,
	Array2Xi&& edgeVertexIndices,
	ArrayX<char>&& nDomainCellEdges,
	Array4Xi&& domainCellEdges,
	const Ref<const Array4Xi>& domainCellVertices
) :
	Mesh ( edgesAndCellsToMesh(
		std::move(vertexCoords),
		std::move(edgeVertexIndices),
		std::move(nDomainCellEdges),
		std::move( domainCellEdges),
		domainCellVertices
	) )
{}


/**
 * @brief Constructor for unstructured meshes. All other constructors chain to 
 * this one, except for the one for full import and the move constructor.
 * Use this ctor when you have
 * vertex coordinates,
 * vertex indices of edges,
 * indices of cells adjacent to each edge,
 * domain cells defined by edge indices (which also requires the number of edges
 * for each cell) and
 * vertex indices for cells. The latter are required to determine cell areas and
 * centroids.
 * Cell centroids are also required to determine the orientation of edge normal
 * vectors.
 * 
 * @param vertexCoords 
 * @param nInnerEdges 
 * @param edgeVertexIndices 
 * @param adjacentCells 
 * @param nDomainCellEdges 
 * @param domainCellEdges 
 * @param domainCellVertices 
 */
Mesh::Mesh(
	Matrix2Xs&& vertexCoords,
	Index nInnerEdges,
	Array2Xi&& edgeVertexIndices,
	Array2Xi&& adjacentCells,
	ArrayX<char>&& nDomainCellEdges,
	Array4Xi&& domainCellEdges,
	const Ref<const Array4Xi>& domainCellVertices
) :
	m_vertices { std::move(vertexCoords) }
{
	m_edges.storeData(
		nInnerEdges,
		std::move(edgeVertexIndices),
		std::move(adjacentCells)
	);
	m_cells.storeData(
		std::move(nDomainCellEdges),
		std::move(domainCellEdges),
		domainCellVertices
	);
	m_edges.storeNormals();

}


auto Mesh::vertices() const -> const Matrix2Xs& {
	return m_vertices;
}


auto Mesh::vertexImp(Index i) const -> Matrix2Xs::ConstColXpr {
	return m_vertices.col(i);
}

/**
 * @brief Creates and returns an object representing an edge, through which the
 * properties of that edge can be easily accessed. The object itself doesn't 
 * hold any data apart from the index passed in and a reference to the edges.
 * @see Mesh::cell()
 * @param i the edge index
 */
auto Mesh::edgeImp(Index i) const -> Edge {
	return Edge { m_edges, i };
}

/**
 * @brief Creates and returns an object representing a cell, through which the
 * properties of that cell can be accessed. The object itself is only a thin
 * wrapper and the performance impact for -O2 or higher is tested to be zero.
 * @see Mesh::edge()
 * @param i the cell index
 */
auto Mesh::cellImp(Index i) const -> Cell {
	return Cell { m_cells, i };
}


auto Mesh::nVerticesImp     () const -> Index {
	return m_vertices.cols();
}
auto Mesh::nEdgesImp        () const -> Index {
	return m_edges.nEdges();
}
auto Mesh::nInnerEdgesImp   () const -> Index {
	return m_edges.nInnerEdges();
}
auto Mesh::nBoundaryEdgesImp() const -> Index {
	return m_edges.nBoundaryEdges();
}
auto Mesh::nCellsImp        () const -> Index {
	return m_cells.nCells();
}
auto Mesh::nDomainCellsImp  () const -> Index {
	return m_cells.nDomainCells();
}
auto Mesh::nGhostCellsImp   () const -> Index {
	return m_cells.nGhostCells();
}

auto Mesh::domainAreaImp    () const -> scalar {
	scalar sum {0};
	OMP_PRAGMA(parallel reduction (+:sum))
	{
		auto [colBeg, nCols] = ompSegment( nDomainCells() );
		sum += m_cells.areas().segment(colBeg, nCols).sum();
	}
	return sum;
}

void Mesh::checkBoundaryPatch( const ArrayXi& ) const {}
// void Mesh::checkBoundaryPatch( const ArrayXi& edges ) const {
// 	Index firstCell { this->edge(edges[0]).ghostCell() };
// 	/* the ghost cells must be consecutively numbered */
// 	for ( Index i{1}; i<edges.size(); ++i ){
// 		if ( firstCell + i != edge(edges[i]).ghostCell() )
// 			throw MeshError(
// 				"Unstructured mesh: Boundary patch must be continuous!"
// 			);
// 	}
// }


} // end namespace hms